#undef DEBUG
module cldwat
!----------------------------------------------------------------------- 
! 
! Purpose: Prognostic cloud water data and methods.
! 
! Public interfaces:
!
! inimc -- Initialize constants
! pcond -- Calculate prognostic condensate
!
! Author: P. Rasch, with Modifications by Minghua Zhang
! January 2010, modified by J. Kay to add precip fluxes for COSP simulator
! 
!-----------------------------------------------------------------------
   use shr_kind_mod,   only: r8 => shr_kind_r8
   use spmd_utils,     only: masterproc
   use ppgrid,         only: pcols, pver, pverp
   use physconst,      only: latvap, latice, cpair
   use cam_abortutils, only: endrun
   use cam_logfile,    only: iulog
   use ref_pres,       only: top_lev => trop_cloud_top_lev

   implicit none

!-----------------------------------------------------------------------
! PUBLIC: Make default data and interfaces private
!-----------------------------------------------------------------------
   private
   save
   public inimc, pcond          ! Public interfaces
   public :: cldwat_init
   integer, public::  ktop      ! Level above 10 hPa

   real(r8),public ::  icritc               ! threshold for autoconversion of cold ice
   real(r8),public ::  icritw               ! threshold for autoconversion of warm ice
!!$   real(r8),public,parameter::  conke  = 1.e-6    ! tunable constant for evaporation of precip
!!$   real(r8),public,parameter::  conke  =  2.e-6    ! tunable constant for evaporation of precip
   real(r8),public ::  conke                ! tunable constant for evaporation of precip
   real(r8),public ::  r3lcrit              ! critical radius where liq conversion begins

!-----------------------------------------------------------------------
! PRIVATE: Everything else is private to this module
!-----------------------------------------------------------------------
   real(r8), private:: rhonot   ! air density at surface
   real(r8), private:: t0       ! Freezing temperature
   real(r8), private:: cldmin   ! assumed minimum cloud amount
   real(r8), private:: small    ! small number compared to unity
   real(r8), private:: c        ! constant for graupel like snow cm**(1-d)/s
   real(r8), private:: d        ! constant for graupel like snow
   real(r8), private:: esi      ! collection efficient for ice by snow
   real(r8), private:: esw      ! collection efficient for water by snow
   real(r8), private:: nos      ! particles snow / cm**4
   real(r8), private:: pi       ! Mathematical constant
   real(r8), private:: gravit   ! Gravitational acceleration at surface
   real(r8), private:: rh2o
   real(r8), private:: prhonos
   real(r8), private:: thrpd    ! numerical three added to d
   real(r8), private:: gam3pd   ! gamma function on (3+d)
   real(r8), private:: gam4pd   ! gamma function on (4+d)
   real(r8), private:: rhoi     ! ice density
   real(r8), private:: rhos     ! snow density
   real(r8), private:: rhow     ! water density
   real(r8), private:: mcon01   ! constants used in cloud microphysics
   real(r8), private:: mcon02   ! constants used in cloud microphysics
   real(r8), private:: mcon03   ! constants used in cloud microphysics
   real(r8), private:: mcon04   ! constants used in cloud microphysics
   real(r8), private:: mcon05   ! constants used in cloud microphysics
   real(r8), private:: mcon06   ! constants used in cloud microphysics
   real(r8), private:: mcon07   ! constants used in cloud microphysics
   real(r8), private:: mcon08   ! constants used in cloud microphysics


! Parameters used in findmcnew
   real(r8) :: capnsi               ! sea ice cloud particles / cm3
   real(r8) :: capnc                ! cold and oceanic cloud particles / cm3
   real(r8) :: capnw                ! warm continental cloud particles / cm3
   real(r8) :: kconst               ! const for terminal velocity (stokes regime)
   real(r8) :: effc                 ! collection efficiency
   real(r8) :: alpha                ! ratio of 3rd moment radius to 2nd
   real(r8) :: capc                 ! constant for autoconversion
   real(r8) :: convfw               ! constant used for fall velocity calculation
   real(r8) :: cracw                ! constant used for rain accreting water
   real(r8) :: critpr               ! critical precip rate collection efficiency changes
   real(r8) :: ciautb               ! coefficient of autoconversion of ice (1/s)
   real(r8) :: psrhmin              ! condensation threshold in polar stratosphere
   logical  :: do_psrhmin

#ifdef DEBUG
   integer, private,parameter ::  nlook = 1  ! Number of points to examine
   integer, private ::  ilook(nlook)         ! Longitude index to examine
   integer, private ::  latlook(nlook)       ! Latitude index to examine
   integer, private ::  lchnklook(nlook)     ! Chunk index to examine
   integer, private ::  icollook(nlook)      ! Column index to examine
#endif

contains
!===============================================================================
subroutine cldwat_init(icritw_in, icritc_in, conke_in, r3lcrit_in, psrhmin_in, do_psrhmin_in )

    real(r8), intent(in) :: icritw_in    !   icritw  = threshold for autoconversion of warm ice  
    real(r8), intent(in) :: icritc_in    !   icritc  = threshold for autoconversion of cold ice  
    real(r8), intent(in) :: conke_in     !   conke   = tunable constant for evaporation of precip
    real(r8), intent(in) :: r3lcrit_in   !   r3lcrit = critical radius where liq conversion begins
    real(r8), intent(in) :: psrhmin_in   ! condensation threadhold in polar stratosphere
    logical,  intent(in) :: do_psrhmin_in

    icritw  = icritw_in
    icritc  = icritc_in
    conke   = conke_in
    r3lcrit = r3lcrit_in
    psrhmin = psrhmin_in
    do_psrhmin = do_psrhmin_in

  end subroutine cldwat_init

subroutine inimc( tmeltx, rhonotx, gravitx, rh2ox)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! initialize constants for the prognostic condensate
! 
! Author: P. Rasch, April 1997
! 
!-----------------------------------------------------------------------
   use pmgrid,       only: plev, plevp
   use ref_pres,     only: pref_mid

   integer k
   real(r8), intent(in) :: tmeltx
   real(r8), intent(in) :: rhonotx
   real(r8), intent(in) :: gravitx
   real(r8), intent(in) :: rh2ox

#ifdef UNICOSMP
   real(r8) signgam              ! variable required by cray gamma function
   external gamma
#endif

   rhonot = rhonotx             ! air density at surface (gm/cm3)
   gravit = gravitx
   rh2o   = rh2ox
   rhos = .1_r8                 ! assumed snow density (gm/cm3)
   rhow = 1._r8                 ! water density
   rhoi = 1._r8                 ! ice density
   esi = 1.0_r8                 ! collection efficient for ice by snow
   esw = 0.1_r8                 ! collection efficient for water by snow
   t0 = tmeltx                  ! approximate freezing temp
   cldmin = 0.02_r8             ! assumed minimum cloud amount
   small = 1.e-22_r8            ! a small number compared to unity
   c = 152.93_r8                ! constant for graupel like snow cm**(1-d)/s
   d = 0.25_r8                  ! constant for graupel like snow
   nos = 3.e-2_r8               ! particles snow / cm**4
   pi = 4._r8*atan(1.0_r8)
   prhonos = pi*rhos*nos
   thrpd = 3._r8 + d
   if (d==0.25_r8) then
      gam3pd = 2.549256966718531_r8 ! only right for d = 0.25
      gam4pd = 8.285085141835282_r8
   else
#ifdef UNICOSMP
      call gamma(3._r8+d, signgam, gam3pd)
      gam3pd = sign(exp(gam3pd),signgam)
      call gamma(4._r8+d, signgam, gam4pd)
      gam4pd = sign(exp(gam4pd),signgam)
      write(iulog,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd
#else
      call endrun(' can only use d ne 0.25 on a cray ')
#endif
   endif
   mcon01 = pi*nos*c*gam3pd/4._r8
   mcon02 = 1._r8/(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4._r8)))
   mcon03 = -(0.5_r8+d/4._r8)
   mcon04 = 4._r8/(4._r8+d)
   mcon05 = (3+d)/(4+d)
   mcon06 = (3+d)/4._r8
   mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06
   mcon08 = -0.5_r8/(4._r8+d)

   if( masterproc ) write(iulog,*) 'cloud water initialization by inimc complete '

! Initialize parameters used by findmcnew
   capnw = 400._r8              ! warm continental cloud particles / cm3
   capnc = 150._r8              ! cold and oceanic cloud particles / cm3
!  capnsi = 40._r8              ! sea ice cloud particles density  / cm3
   capnsi = 75._r8              ! sea ice cloud particles density  / cm3

   kconst = 1.18e6_r8           ! const for terminal velocity

!  effc = 1._r8                 ! autoconv collection efficiency following boucher 96
!  effc = .55*0.05_r8           ! autoconv collection efficiency following baker 93
   effc = 0.55_r8               ! autoconv collection efficiency following tripoli and cotton
!  effc = 0._r8    ! turn off warm-cloud autoconv
   alpha = 1.1_r8**4
   capc = pi**(-.333_r8)*kconst*effc *(0.75_r8)**(1.333_r8)*alpha  ! constant for autoconversion

! critical precip rate at which we assume the collector drops can change the
! drop size enough to enhance the auto-conversion process (mm/day)
   critpr = 0.5_r8
   convfw = 1.94_r8*2.13_r8*sqrt(rhow*1000._r8*9.81_r8*2.7e-4_r8)

! liquid microphysics
!  cracw = 6_r8                 ! beheng
   cracw = .884_r8*sqrt(9.81_r8/(rhow*1000._r8*2.7e-4_r8)) ! tripoli and cotton

! ice microphysics
   ciautb = 5.e-4_r8

   if ( masterproc ) then
      write(iulog,*)'tuning parameters cldwat: icritw',icritw,'icritc',icritc,'conke',conke,'r3lcrit',r3lcrit
      write(iulog,*)'tuning parameters cldwat: capnw',capnw,'capnc',capnc,'capnsi',capnsi,'kconst',kconst
      write(iulog,*)'tuning parameters cldwat: effc',effc,'alpha',alpha,'capc',capc
      write(iulog,*)'tuning parameters cldwat: critpr',critpr,'convfw',convfw,'cracw',cracw,'ciautb',ciautb
   endif

end subroutine inimc

subroutine pcond (lchnk   ,ncol    ,troplev ,dlat    , &
                  tn      ,ttend   ,qn      ,qtend   ,omega   , &
                  cwat    ,p       ,pdel    ,cldn    ,fice    , fsnow, &
                  cme     ,prodprec,prodsnow,evapprec,evapsnow,evapheat, prfzheat, &     
                  meltheat,precip  ,snowab  ,deltat  ,fwaut   , &
                  fsaut   ,fracw   ,fsacw   ,fsaci   ,lctend  , &
                  rhdfda  ,rhu00   ,landm   ,seaicef ,zi      ,ice2pr, liq2pr, &
                  liq2snow, snowh, rkflxprc, rkflxsnw, pracwo, psacwo, psacio)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! The public interface to the cloud water parameterization
! returns tendencies to water vapor, temperature and cloud water variables
! 
! For basic method 
!  See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
!  model climate using diagnosed and 
!  predicted condensate parameterizations, 1998, J. Clim., 11,
!  pp1587---1614.
! 
! For important modifications to improve the method of determining
! condensation/evaporation see Zhang et al (2001, in preparation)
!
! Authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson
!          B. A. Boville (latent heat of fusion)
!-----------------------------------------------------------------------
   use wv_saturation, only: qsat, estblf, svp_to_qsat, findsp_vc
   use physconst, only: epsilo
!
!---------------------------------------------------------------------
!
! Input Arguments
!
   integer,  intent(in) :: lchnk                ! chunk identifier
   integer,  intent(in) :: ncol                 ! number of atmospheric columns
   integer,  intent(in) :: troplev(pcols)       ! tropopause level
   real(r8), intent(in) :: dlat(pcols)          ! latitudes in degrees
   real(r8), intent(in) :: fice(pcols,pver)     ! fraction of cwat that is ice
   real(r8), intent(in) :: fsnow(pcols,pver)    ! fraction of rain that freezes to snow
   real(r8), intent(in) :: cldn(pcols,pver)     ! new value of cloud fraction    (fraction)
   real(r8), intent(in) :: cwat(pcols,pver)     ! cloud water (kg/kg)
   real(r8), intent(in) :: omega(pcols,pver)    ! vert pressure vel (Pa/s)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure          (K)
   real(r8), intent(in) :: pdel(pcols,pver)     ! pressure thickness (Pa)
   real(r8), intent(in) :: qn(pcols,pver)       ! new water vapor    (kg/kg)
   real(r8), intent(in) :: qtend(pcols,pver)    ! mixing ratio tend  (kg/kg/s)
   real(r8), intent(in) :: tn(pcols,pver)       ! new temperature    (K)
   real(r8), intent(in) :: ttend(pcols,pver)    ! temp tendencies    (K/s)
   real(r8), intent(in) :: deltat               ! time step to advance solution over
   real(r8), intent(in) :: lctend(pcols,pver)   ! cloud liquid water tendencies   ====wlin
   real(r8), intent(in) :: rhdfda(pcols,pver)   ! dG(a)/da, rh=G(a), when rh>u00  ====wlin
   real(r8), intent(in) :: rhu00 (pcols,pver)   ! Rhlim for cloud                 ====wlin
   real(r8), intent(in) :: landm(pcols)         ! Land fraction ramped over water (fraction)
   real(r8), intent(in) :: seaicef(pcols)       ! sea ice fraction  (fraction)
   real(r8), intent(in) :: zi(pcols,pverp)      ! layer interfaces (m)
   real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
!
! Output Arguments
!
   real(r8), intent(out) :: cme     (pcols,pver) ! rate of cond-evap of condensate (1/s)
   real(r8), intent(out) :: prodprec(pcols,pver) ! rate of conversion of condensate to precip (1/s)
   real(r8), intent(out) :: evapprec(pcols,pver) ! rate of evaporation of falling precip (1/s)
   real(r8), intent(out) :: evapsnow(pcols,pver) ! rate of evaporation of falling snow (1/s)
   real(r8), intent(out) :: evapheat(pcols,pver) ! heating rate due to evaporation of precip (W/kg)
   real(r8), intent(out) :: prfzheat(pcols,pver) ! heating rate due to freezing of precip (W/kg)
   real(r8), intent(out) :: meltheat(pcols,pver) ! heating rate due to snow melt (W/kg)
   real(r8), intent(out) :: precip(pcols)        ! rate of precipitation (kg / (m**2 * s))
   real(r8), intent(out) :: snowab(pcols)        ! rate of snow (kg / (m**2 * s))
   real(r8), intent(out) :: ice2pr(pcols,pver)   ! rate of conversion of ice to precip
   real(r8), intent(out) :: liq2pr(pcols,pver)   ! rate of conversion of liquid to precip
   real(r8), intent(out) :: liq2snow(pcols,pver) ! rate of conversion of liquid to snow
   real(r8), intent(out) :: rkflxprc(pcols,pverp)   ! grid-box mean RK flux_large_scale_cloud_rain+snow at interfaces (kg m^-2 s^-1)
   real(r8), intent(out) :: rkflxsnw(pcols,pverp)   ! grid-box mean RK flux_large_scale_cloud_snow at interfaces (kg m^-2 s^-1)
! intent(out)s here for pcond to pass to stratiform.F90 to be addflded/outflded
   real(r8), intent(out) :: pracwo(pcols,pver)      ! accretion of cloud water by rain (1/s)
   real(r8), intent(out) :: psacwo(pcols,pver)      ! accretion of cloud water by snow (1/s)
   real(r8), intent(out) :: psacio(pcols,pver)      ! accretion of cloud ice by snow (1/s)

   real(r8) nice2pr     ! rate of conversion of ice to snow
   real(r8) nliq2pr     ! rate of conversion of liquid to precip
   real(r8) nliq2snow   ! rate of conversion of liquid to snow
   real(r8) prodsnow(pcols,pver) ! rate of production of snow

!
! Local workspace
!
   real(r8) :: precab(pcols)        ! rate of precipitation (kg / (m**2 * s))
   integer i                 ! work variable
   integer iter              ! #iterations for precipitation calculation
   integer k                 ! work variable
   integer l                 ! work variable

   real(r8) cldm(pcols)          ! mean cloud fraction over the time step
   real(r8) cldmax(pcols)        ! max cloud fraction above
   real(r8) coef(pcols)          ! conversion time scale for condensate to rain
   real(r8) cwm(pcols)           ! cwat mixing ratio at midpoint of time step
   real(r8) cwn(pcols)           ! cwat mixing ratio at end
   real(r8) denom                ! work variable
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) es(pcols)            ! sat. vapor pressure
   real(r8) fracw(pcols,pver)    ! relative importance of collection of liquid by rain
   real(r8) fsaci(pcols,pver)    ! relative importance of collection of ice by snow
   real(r8) fsacw(pcols,pver)    ! relative importance of collection of liquid by snow
   real(r8) fsaut(pcols,pver)    ! relative importance of ice auto conversion
   real(r8) fwaut(pcols,pver)    ! relative importance of warm cloud autoconversion
   real(r8) gamma(pcols)         ! d qs / dT
   real(r8) icwc(pcols)          ! in-cloud water content (kg/kg)
   real(r8) mincld               ! a small cloud fraction to avoid / zero
   real(r8),parameter ::omsm=0.99999_r8                 ! a number just less than unity (for rounding)
   real(r8) prprov(pcols)        ! provisional value of precip at btm of layer
   real(r8) prtmp                ! work variable
   real(r8) q(pcols,pver)        ! mixing ratio before time step ignoring condensate
   real(r8) qs(pcols)            ! spec. hum. of water vapor
   real(r8) qsn, esn             ! work variable
   real(r8) qsp(pcols,pver)      ! sat pt mixing ratio
   real(r8) qtl(pcols)           ! tendency which would saturate the grid box in deltat
   real(r8) qtmp, ttmp           ! work variable
   real(r8) relhum1(pcols)        ! relative humidity
   real(r8) relhum(pcols)        ! relative humidity
!!$   real(r8) tc                   ! crit temp of transition to ice
   real(r8) t(pcols,pver)        ! temp before time step ignoring condensate
   real(r8) tsp(pcols,pver)      ! sat pt temperature
   real(r8) pol                  ! work variable
   real(r8) cdt                  ! work variable
   real(r8) wtthick              ! work variable

! Extra local work space for cloud scheme modification       

   real(r8) cpohl                !cpair/Latvap
   real(r8) hlocp                !Latvap/cpair
   real(r8) dto2                 !0.5*deltat (delta=2.0*dt)
   real(r8) calpha(pcols)        !alpha of new C - E scheme formulation
   real(r8) cbeta (pcols)        !beta  of new C - E scheme formulation
   real(r8) cbetah(pcols)        !beta_hat at saturation portion 
   real(r8) cgamma(pcols)        !gamma of new C - E scheme formulation
   real(r8) cgamah(pcols)        !gamma_hat at saturation portion
   real(r8) rcgama(pcols)        !gamma/gamma_hat
   real(r8) csigma(pcols)        !sigma of new C - E scheme formulation
   real(r8) cmec1 (pcols)        !c1    of new C - E scheme formulation
   real(r8) cmec2 (pcols)        !c2    of new C - E scheme formulation
   real(r8) cmec3 (pcols)        !c3    of new C - E scheme formulation
   real(r8) cmec4 (pcols)        !c4    of new C - E scheme formulation
   real(r8) cmeres(pcols)        !residual cond of over-sat after cme and evapprec
   real(r8) ctmp                 !a scalar representation of cmeres
   real(r8) clrh2o               ! Ratio of latvap to water vapor gas const
   real(r8) ice(pcols,pver)    ! ice mixing ratio
   real(r8) liq(pcols,pver)    ! liquid mixing ratio
   real(r8) rcwn(pcols,2,pver), rliq(pcols,2,pver), rice(pcols,2,pver)
   real(r8) cwnsave(pcols,2,pver), cmesave(pcols,2,pver)
   real(r8) prodprecsave(pcols,2,pver)
   logical error_found

   real(r8) :: rhu_adj(pcols,pver)   ! adjusted rhlim for dehydration 
!
!------------------------------------------------------------
!
   clrh2o = latvap/rh2o   ! Ratio of latvap to water vapor gas const
#ifdef PERGRO
   mincld = 1.e-4_r8
   iter = 1   ! number of times to iterate the precipitation calculation
#else
   mincld = 1.e-4_r8
   iter = 2
#endif
!   omsm = 0.99999
   cpohl = cpair/latvap
   hlocp = latvap/cpair
   dto2=0.5_r8*deltat
!
! Constant for computing rate of evaporation of precipitation:
!
!!$   conke = 1.e-5
!!$   conke = 1.e-6
!
! initialize a few single level fields
!
   do i = 1,ncol
      precip(i) = 0.0_r8
      precab(i) = 0.0_r8
      snowab(i) = 0.0_r8
      cldmax(i) = 0.0_r8
   end do
!
! initialize multi-level fields 
!
   do k = 1,pver
      do i = 1,ncol
         q(i,k) = qn(i,k) 
         t(i,k) = tn(i,k)
!         q(i,k)=qn(i,k)-qtend(i,k)*deltat
!         t(i,k)=tn(i,k)-ttend(i,k)*deltat
    end do
   end do
   cme     (:ncol,:) = 0._r8
   evapprec(:ncol,:) = 0._r8
   prodprec(:ncol,:) = 0._r8
   evapsnow(:ncol,:) = 0._r8
   prodsnow(:ncol,:) = 0._r8
   evapheat(:ncol,:) = 0._r8
   meltheat(:ncol,:) = 0._r8
   prfzheat(:ncol,:) = 0._r8
   ice2pr(:ncol,:)   = 0._r8
   liq2pr(:ncol,:)   = 0._r8
   liq2snow(:ncol,:) = 0._r8
   fwaut(:ncol,:) = 0._r8
   fsaut(:ncol,:) = 0._r8
   fracw(:ncol,:) = 0._r8
   fsacw(:ncol,:) = 0._r8
   fsaci(:ncol,:) = 0._r8
   rkflxprc(:ncol,:) = 0._r8
   rkflxsnw(:ncol,:) = 0._r8

   pracwo(:ncol,:) = 0._r8
   psacwo(:ncol,:) = 0._r8
   psacio(:ncol,:) = 0._r8
!
! find the wet bulb temp and saturation value
! for the provisional t and q without condensation
!
   do 800 k = top_lev,pver

      ! "True" means that ice will be taken into account.
      call findsp_vc(qn(:ncol,k), tn(:ncol,k), p(:ncol,k), .true., &
           tsp(:ncol,k), qsp(:ncol,k))

      call qsat(t(:ncol,k), p(:ncol,k), &
           es(:ncol), qs(:ncol), gam=gamma(:ncol))
      do i = 1,ncol
         relhum(i) = q(i,k)/qs(i)
!
         cldm(i) = max(cldn(i,k),mincld)
!
! the max cloud fraction above this level
!
         cldmax(i) = max(cldmax(i), cldm(i))

! define the coefficients for C - E calculation

         calpha(i) = 1.0_r8/qs(i)
         cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl
         cbetah(i) = 1.0_r8/qs(i)*gamma(i)*cpohl
         cgamma(i) = calpha(i)+latvap*cbeta(i)/cpair
         cgamah(i) = calpha(i)+latvap*cbetah(i)/cpair
         rcgama(i) = cgamma(i)/cgamah(i)

         if(cldm(i) > mincld) then
            icwc(i) = max(0._r8,cwat(i,k)/cldm(i))
         else
            icwc(i) = 0.0_r8
         endif
!PJR the above logic give zero icwc with nonzero cwat, dont like it!
!PJR generates problems with csigma
!PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
!         icwc(i) = max(1.e-8_r8,cwat(i,k)/cldm(i))

!
! initial guess of evaporation, will be updated within iteration
!
         evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
                        *(1._r8 - min(relhum(i),1._r8))

!
! zero cmeres before iteration for each level
!
         cmeres(i)=0.0_r8

      end do
      do i = 1,ncol
!
! fractions of ice at this level
!
!!$         tc = t(i,k) - t0
!!$         fice(i,k) = max(0._r8,min(-tc*0.05,1.0_r8))
!
! calculate the cooling due to a phase change of the rainwater
! from above
!
         if (t(i,k) >= t0) then
            meltheat(i,k) =  -latice * snowab(i) * gravit/pdel(i,k)
            snowab(i) = 0._r8
         else
            meltheat(i,k) = 0._r8
         endif
      end do

!
! calculate cme and formation of precip. 
!
! The cloud microphysics is highly nonlinear and coupled with cme
! Both rain processes and cme are calculated iteratively.
! 
      do 100 l = 1,iter

        do i = 1,ncol

!
! calculation of cme has 4 scenarios
! ==================================
!
           call relhum_min_adj( ncol, troplev, dlat, rhu00, rhu_adj )

           if(relhum(i) > rhu_adj(i,k)) then
    
           ! 1. whole grid saturation
           ! ========================
              if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then
                 cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i)

           ! 2. fractional saturation
           ! ========================
              else
                 if (rhdfda(i,k) .eq. 0._r8 .and. icwc(i).eq.0._r8) then
                    write (iulog,*) ' cldwat.F90:  empty rh cloud ', i, k, lchnk
                    write (iulog,*) ' relhum, iter ', relhum(i), l, rhu_adj(i,k), cldm(i), cldn(i,k)
                    call endrun ()
                 endif
                  csigma(i) = 1.0_r8/(rhdfda(i,k)+cgamma(i)*icwc(i))
                  cmec1(i) = (1.0_r8-cldm(i))*csigma(i)*rhdfda(i,k)
                  cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0_r8-rcgama(i)*cldm(i))*   &
                             csigma(i)*calpha(i)*icwc(i)
                  cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) +  &
                           (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i)
                  cmec4(i) = csigma(i)*cgamma(i)*icwc(i)

                  ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er
  
                  cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k)  &
                             -cmec3(i)*ttend(i,k) + cmec4(i)*evapprec(i,k)
               endif

           ! 3. when rh < rhu00, evaporate existing cloud water
           ! ================================================== 
           else if(cwat(i,k) > 0.0_r8)then
              ! liquid water should be evaporated but not to exceed 
              ! saturation point. if qn > qsp, not to evaporate cwat
              cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat 

           ! 4. no condensation nor evaporation
           ! ==================================
           else
              cme(i,k)=0.0_r8
           endif

  
        end do    !end loop for cme update

! Because of the finite time step, 
! place a bound here not to exceed wet bulb point
! and not to evaporate more than available water
!
         do i = 1, ncol
            qtmp = qn(i,k) - cme(i,k)*deltat

! possibilities to have qtmp > qsp
!
!   1. if qn > qs(tn), it condenses; 
!      if after applying cme,  qtmp > qsp,  more condensation is applied. 
!      
!   2. if qn < qs, evaporation should not exceed qsp,
    
            if(qtmp > qsp(i,k)) then
              cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat
            endif

!
! if net evaporation, it should not exceed available cwat
!
            if(cme(i,k) < -cwat(i,k)/deltat)  &
               cme(i,k) = -cwat(i,k)/deltat
!
! addition of residual condensation from previous step of iteration
!
            cme(i,k) = cme(i,k) + cmeres(i)

         end do

         !      limit cme for roundoff errors
         do i = 1, ncol
            cme(i,k) = cme(i,k)*omsm
         end do

         do i = 1,ncol
!
! as a safe limit, condensation should not reduce grid mean rh below rhu00
! 
           if(cme(i,k) > 0.0_r8 .and. relhum(i) > rhu_adj(i,k) )  &
              cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu_adj(i,k))/deltat)
!
! initial guess for cwm (mean cloud water over time step) if 1st iteration
!
           if(l < 2) then
             cwm(i) = max(cwat(i,k)+cme(i,k)*dto2,  0._r8)
           endif

         enddo

! provisional precipitation falling through model layer
         do i = 1,ncol
!!$            prprov(i) =  precab(i) + prodprec(i,k)*pdel(i,k)/gravit
! rain produced in this layer not too effective in collection process
            wtthick = max(0._r8,min(0.5_r8,((zi(i,k)-zi(i,k+1))/1000._r8)**2))
            prprov(i) =  precab(i) + wtthick*prodprec(i,k)*pdel(i,k)/gravit
         end do

! calculate conversion of condensate to precipitation by cloud microphysics 
         call findmcnew (lchnk   ,ncol    , &
                         k       ,prprov  ,snowab,  t       ,p        , &
                         cwm     ,cldm    ,cldmax  ,fice(1,k),coef    , &
                         fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), &
                         landm, seaicef, snowh, pracwo(1,k), psacwo(1,k), psacio(1,k))

!
! calculate the precip rate
!
         error_found = .false.
         do i = 1,ncol
            if (cldm(i) > 0) then  
!
! first predict the cloud water
!
               cdt = coef(i)*deltat
               if (cdt > 0.01_r8) then
                  pol = cme(i,k)/coef(i) ! production over loss
                  cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol)
               else
                  cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt))
               endif
!
! now back out the tendency of net rain production
!
               prodprec(i,k) =  max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat)
            else
               prodprec(i,k) = 0.0_r8
               cwn(i) = 0._r8
            endif

            ! provisional calculation of conversion terms
            ice2pr(i,k) = prodprec(i,k)*(fsaut(i,k)+fsaci(i,k))
            liq2pr(i,k) = prodprec(i,k)*(fwaut(i,k)+fsacw(i,k)+fracw(i,k))
!old        liq2snow(i,k) = prodprec(i,k)*fsacw(i,k)

!           revision suggested by Jim McCaa
!           it controls the amount of snow hitting the sfc 
!           by forcing a lot of conversion of cloud liquid to snow phase
!           it might be better done later by an explicit representation of 
!           rain accreting ice (and freezing), or by an explicit freezing of raindrops
            liq2snow(i,k) = max(prodprec(i,k)*fsacw(i,k), fsnow(i,k)*liq2pr(i,k))

            ! bounds
            nice2pr = min(ice2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*fice(i,k)/deltat)
            nliq2pr = min(liq2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*(1._r8-fice(i,k))/deltat)
!            write(iulog,*) ' prodprec ', i, k, prodprec(i,k)
!            write(iulog,*) ' nliq2pr, nice2pr ', nliq2pr, nice2pr
            if (liq2pr(i,k).ne.0._r8) then
               nliq2snow = liq2snow(i,k)*nliq2pr/liq2pr(i,k)   ! correction
            else
               nliq2snow = liq2snow(i,k)
            endif

!           avoid roundoff problems generating negatives
            nliq2snow = nliq2snow*omsm
            nliq2pr = nliq2pr*omsm
            nice2pr = nice2pr*omsm
            
!           final estimates of conversion to precip and snow
            prodprec(i,k) = (nliq2pr + nice2pr)
            prodsnow(i,k) = (nice2pr + nliq2snow)

            rcwn(i,l,k) =  cwat(i,k) + (cme(i,k)-   prodprec(i,k))*deltat
            rliq(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)*(1._r8-fice(i,k)) - nliq2pr * deltat
            rice(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)* fice(i,k)      -    nice2pr                     *deltat

!           Save for sanity check later...  
!           Putting sanity checks inside loops 100 and 800 screws up the 
!           IBM compiler for reasons as yet unknown.  TBH
            cwnsave(i,l,k)      = cwn(i)
            cmesave(i,l,k)      = cme(i,k)
            prodprecsave(i,l,k) = prodprec(i,k)
!           End of save for sanity check later...  

!           final version of condensate to precip terms
            liq2pr(i,k) = nliq2pr
            liq2snow(i,k) = nliq2snow
            ice2pr(i,k) = nice2pr

            cwn(i) = rcwn(i,l,k)
!
! update any remaining  provisional values
!
            cwm(i) = (cwn(i) + cwat(i,k))*0.5_r8
!
! update in cloud water
!
            if(cldm(i) > mincld) then
               icwc(i) = cwm(i)/cldm(i)
            else
               icwc(i) = 0.0_r8
            endif
!PJR the above logic give zero icwc with nonzero cwat, dont like it!
!PJR generates problems with csigma
!PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
!         icwc(i) = max(1.e-8_r8,cwm(i)/cldm(i))

         end do              ! end of do i = 1,ncol

!
! calculate provisional value of cloud water for
! evaporation of precipitate (evapprec) calculation
!
      do i = 1,ncol
         qtmp = qn(i,k) - cme(i,k)*deltat
         ttmp = tn(i,k) + deltat/cpair * ( meltheat(i,k)       &
              + (latvap + latice*fice(i,k)) * cme(i,k) )
         esn = estblf(ttmp)
         qsn = svp_to_qsat(esn, p(i,k))
         qtl(i) = max((qsn - qtmp)/deltat,0._r8)
         relhum1(i) = qtmp/qsn
      end do
!
      do i = 1,ncol
#ifdef PERGRO
         evapprec(i,k) = conke*(1._r8 - max(cldm(i),mincld))* &
                         sqrt(precab(i))*(1._r8 - min(relhum1(i),1._r8))
#else
         evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
                         *(1._r8 - min(relhum1(i),1._r8))
#endif
!
! limit the evaporation to the amount which is entering the box
! or saturates the box
!
         prtmp = precab(i)*gravit/pdel(i,k)
         evapprec(i,k) = min(evapprec(i,k), prtmp, qtl(i))*omsm
#ifdef PERGRO
!           zeroing needed for pert growth
         evapprec(i,k) = 0._r8
#endif
!
! Partition evaporation of precipitate between rain and snow using
! the fraction of snow falling into the box. Determine the heating
! due to evaporation. Note that evaporation is positive (loss of precip,
! gain of vapor) and that heating is negative.
         if (evapprec(i,k) > 0._r8) then
            evapsnow(i,k) = evapprec(i,k) * snowab(i) / precab(i)
            evapheat(i,k) = -latvap * evapprec(i,k) - latice * evapsnow(i,k)
         else 
            evapsnow(i,k) = 0._r8
            evapheat(i,k) = 0._r8
         end if
! Account for the latent heat of fusion for liquid drops collected by falling snow
         prfzheat(i,k) = latice * liq2snow(i,k)
      end do

! now remove the residual of any over-saturation. Normally,
! the oversaturated water vapor should have been removed by 
! cme formulation plus constraints by wet bulb tsp/qsp
! as computed above. However, because of non-linearity,
! addition of (cme-evapprec) to update t and q may still cause
! a very small amount of over saturation. It is called a
! residual of over-saturation because theoretically, cme
! should have taken care of all of large scale condensation.
! 

       do i = 1,ncol
          qtmp = qn(i,k)-(cme(i,k)-evapprec(i,k))*deltat
          ttmp = tn(i,k) + deltat/cpair * ( meltheat(i,k) + evapheat(i,k) + prfzheat(i,k)      &
              + (latvap + latice*fice(i,k)) * cme(i,k) )

          call qsat(ttmp, p(i,k), esn, qsn, dqsdt=dqsdt)

          if( qtmp > qsn ) then
             !
             !now extra condensation to bring air to just saturation
             !
             ctmp = (qtmp-qsn)/(1._r8+hlocp*dqsdt)/deltat
             cme(i,k) = cme(i,k)+ctmp
!
! save residual on cmeres to addtion to cme on entering next iteration
! cme exit here contain the residual but overrided if back to iteration
!
             cmeres(i) = ctmp
          else
             cmeres(i) = 0.0_r8
          endif
       end do
              
 100 continue              ! end of do l = 1,iter

!
! precipitation
!
      do i = 1,ncol
         precip(i) = precip(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
         precab(i) = precab(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
         if(precab(i).lt.0._r8) precab(i)=0._r8
!         snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodprec(i,k)*fice(i,k) - evapsnow(i,k))
         snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodsnow(i,k) - evapsnow(i,k))

         ! If temperature above freezing, all precip is rain flux.  if temperature below freezing, all precip is snow flux.
         rkflxprc(i,k+1) = precab(i)   !! making this consistent with other precip fluxes.  prc = rain + snow
         !!rkflxprc(i,k+1) = precab(i) - snowab(i)
         rkflxsnw(i,k+1) = snowab(i)

!!$         if ((precab(i)) < 1.e-10) then      
!!$            precab(i) = 0.
!!$            snowab(i) = 0.
!!$         endif
      end do
 800 continue                ! level loop (k=1,pver)

! begin sanity checks
   error_found = .false.
   do k = top_lev,pver
      do l = 1,iter
         do i = 1,ncol
	    if (abs(rcwn(i,l,k)).lt.1.e-300_r8) rcwn(i,l,k) = 0._r8
	    if (abs(rliq(i,l,k)).lt.1.e-300_r8) rliq(i,l,k) = 0._r8
	    if (abs(rice(i,l,k)).lt.1.e-300_r8) rice(i,l,k) = 0._r8
            if (rcwn(i,l,k).lt.0._r8) error_found = .true.
            if (rliq(i,l,k).lt.0._r8) error_found = .true.
            if (rice(i,l,k).lt.0._r8) error_found = .true.
         enddo
      enddo
   enddo
   if (error_found) then
      do k = top_lev,pver
         do l = 1,iter
            do i = 1,ncol
               if (rcwn(i,l,k).lt.0._r8) then
                  write(iulog,*) ' prob with neg rcwn1 ', rcwn(i,l,k),  &
                     cwnsave(i,l,k)
                  write(iulog,*) ' cwat, cme*deltat, prodprec*deltat ', &
                     cwat(i,k), cmesave(i,l,k)*deltat,               &
                     prodprecsave(i,l,k)*deltat,                     &
                     (cmesave(i,l,k)-prodprecsave(i,l,k))*deltat
                  call endrun('PCOND')
               endif
               if (rliq(i,l,k).lt.0._r8) then
                  write(iulog,*) ' prob with neg rliq1 ', rliq(i,l,k)
                  call endrun('PCOND')
               endif
               if (rice(i,l,k).lt.0._r8) then
                  write(iulog,*) ' prob with neg rice ', rice(i,l,k)
                  call endrun('PCOND')
               endif
            enddo
         enddo
      enddo
   end if
! end sanity checks

   return
end subroutine pcond

!##############################################################################

subroutine findmcnew (lchnk   ,ncol    , &
                      k       ,precab  ,snowab,  t       ,p       , &
                      cwm     ,cldm    ,cldmax  ,fice    ,coef    , &
                      fwaut   ,fsaut   ,fracw   ,fsacw   ,fsaci   , &
                      landm   ,seaicef ,snowh   ,pracwo  ,psacwo  ,psacio )
 
!----------------------------------------------------------------------- 
! 
! Purpose: 
! calculate the conversion of condensate to precipitate
! 
! Method: 
! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
!  model climate using diagnosed and 
!  predicted condensate parameterizations, 1998, J. Clim., 11,
!  pp1587---1614.
! 
! Author: P. Rasch
! 
!-----------------------------------------------------------------------
   use phys_grid, only: get_rlat_all_p
!
! input args
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns
   integer, intent(in) :: k                     ! level index

   real(r8), intent(in) :: precab(pcols)        ! rate of precipitation from above (kg / (m**2 * s))
   real(r8), intent(in) :: t(pcols,pver)        ! temperature       (K)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure          (Pa)
   real(r8), intent(in) :: cldm(pcols)          ! cloud fraction
   real(r8), intent(in) :: cldmax(pcols)        ! max cloud fraction above this level
   real(r8), intent(in) :: cwm(pcols)           ! condensate mixing ratio (kg/kg)
   real(r8), intent(in) :: fice(pcols)          ! fraction of cwat that is ice
   real(r8), intent(in) :: landm(pcols)         ! Land fraction ramped over water
   real(r8), intent(in) :: seaicef(pcols)       ! sea ice fraction 
   real(r8), intent(in) :: snowab(pcols)        ! rate of snow from above (kg / (m**2 * s))
   real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)

! output arguments
   real(r8), intent(out) :: coef(pcols)          ! conversion rate (1/s)
   real(r8), intent(out) :: fwaut(pcols)         ! relative importance of liquid autoconversion (a diagnostic)
   real(r8), intent(out) :: fsaut(pcols)         ! relative importance of ice autoconversion    (a diagnostic)
   real(r8), intent(out) :: fracw(pcols)         ! relative importance of rain accreting liquid (a diagnostic)
   real(r8), intent(out) :: fsacw(pcols)         ! relative importance of snow accreting liquid (a diagnostic)
   real(r8), intent(out) :: fsaci(pcols)         ! relative importance of snow accreting ice    (a diagnostic)
   real(r8), intent(out) :: pracwo(pcols)        ! accretion of cloud water by rain (1/s)
   real(r8), intent(out) :: psacwo(pcols)        ! accretion of cloud water by snow (1/s)
   real(r8), intent(out) :: psacio(pcols)        ! accretion of cloud ice by snow (1/s)


! work variables

   integer i
   integer ii
   integer ind(pcols)
   integer ncols

   real(r8), parameter :: degrad = 57.296_r8 ! divide by this to convert degrees to radians
   real(r8) capn                 ! local cloud particles / cm3
   real(r8) capnoice             ! local cloud particles when not over sea ice / cm3
   real(r8) ciaut                ! coefficient of autoconversion of ice (1/s)
   real(r8) cldloc(pcols)        ! non-zero amount of cloud
   real(r8) cldpr(pcols)         ! assumed cloudy volume occupied by rain and cloud
   real(r8) con1                 ! work constant
   real(r8) con2                 ! work constant
   real(r8) csacx                ! constant used for snow accreting liquid or ice
!!$   real(r8) dtice                ! interval for transition from liquid to ice
   real(r8) icemr(pcols)         ! in-cloud ice mixing ratio
   real(r8) icrit                ! threshold for autoconversion of ice
   real(r8) liqmr(pcols)         ! in-cloud liquid water mixing ratio
   real(r8) pracw                ! rate of rain accreting water
   real(r8) prlloc(pcols)        ! local rain flux in mm/day
   real(r8) prscgs(pcols)        ! local snow amount in cgs units
   real(r8) psaci                ! rate of collection of ice by snow (lin et al 1983)
   real(r8) psacw                ! rate of collection of liquid by snow (lin et al 1983)
   real(r8) psaut                ! rate of autoconversion of ice condensate
   real(r8) ptot                 ! total rate of conversion
   real(r8) pwaut                ! rate of autoconversion of liquid condensate
   real(r8) r3l                  ! volume radius
   real(r8) rainmr(pcols)        ! in-cloud rain mixing ratio
   real(r8) rat1                 ! work constant
   real(r8) rat2                 ! work constant
!!$   real(r8) rdtice               ! recipricol of dtice
   real(r8) rho(pcols)           ! density (mks units)
   real(r8) rhocgs               ! density (cgs units)
   real(r8) rlat(pcols)          ! latitude (radians)
   real(r8) snowfr               ! fraction of precipate existing as snow
   real(r8) totmr(pcols)         ! in-cloud total condensate mixing ratio
   real(r8) vfallw               ! fall speed of precipitate as liquid
   real(r8) wp                   ! weight factor used in calculating pressure dep of autoconversion
   real(r8) wsi                  ! weight factor for sea ice
   real(r8) wt                   ! fraction of ice
   real(r8) wland                ! fraction of land

!      real(r8) csaci
!      real(r8) csacw
!      real(r8) cwaut
!      real(r8) efact
!      real(r8) lamdas
!      real(r8) lcrit
!      real(r8) rcwm
!      real(r8) r3lc2
!      real(r8) snowmr(pcols)
!      real(r8) vfalls

   real(r8) ftot

!     inline statement functions
   real(r8) heavy, heavym, a1, a2, heavyp, heavymp
   heavy(a1,a2) = max(0._r8,sign(1._r8,a1-a2))  ! heavyside function
   heavym(a1,a2) = max(0.01_r8,sign(1._r8,a1-a2))  ! modified heavyside function
!
! New heavyside functions to perhaps address error growth problems
!
   heavyp(a1,a2) = a1/(a2+a1+1.e-36_r8)
   heavymp(a1,a2) = (a1+0.01_r8*a2)/(a2+a1+1.e-36_r8)

!
! find all the points where we need to do the microphysics
! and set the output variables to zero
!
   ncols = 0
   do i = 1,ncol
      coef(i) = 0._r8
      fwaut(i) = 0._r8
      fsaut(i) = 0._r8
      fracw(i) = 0._r8
      fsacw(i) = 0._r8
      fsaci(i) = 0._r8
      liqmr(i) = 0._r8
      rainmr(i) = 0._r8
      if (cwm(i) > 1.e-20_r8) then
         ncols = ncols + 1
         ind(ncols) = i
      endif
   end do

   do ii = 1,ncols
      i = ind(ii)
!
! the local cloudiness at this level
!
      cldloc(i) = max(cldmin,cldm(i))
!
! a weighted mean between max cloudiness above, and this layer
!
      cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5_r8)
!
! decompose the suspended condensate into
! an incloud liquid and ice phase component
!
      totmr(i) = cwm(i)/cldloc(i)
      icemr(i) = totmr(i)*fice(i)
      liqmr(i) = totmr(i)*(1._r8-fice(i))
!
! density
!
      rho(i) = p(i,k)/(287._r8*t(i,k))
      rhocgs = rho(i)*1.e-3_r8     ! density in cgs units
!
! decompose the precipitate into a liquid and ice phase
!
      if (t(i,k) > t0) then
         vfallw = convfw/sqrt(rho(i))
         rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i))
         snowfr = 0
!        snowmr(i)
      else
         snowfr = 1
         rainmr(i) = 0._r8
      endif
!     rainmr(i) = (precab(i)-snowab(i))/(rho(i)*vfallw*cldpr(i))
!
! local snow amount in cgs units
!
      prscgs(i) = precab(i)/cldpr(i)*0.1_r8*snowfr
!     prscgs(i) = snowab(i)/cldpr(i)*0.1
!
! local rain amount in mm/day
!
      prlloc(i) = precab(i)*86400._r8/cldpr(i)
   end do

   con1 = 1._r8/(1.333_r8*pi)**0.333_r8 * 0.01_r8 ! meters
!
! calculate the conversion terms
!
   call get_rlat_all_p(lchnk, ncol, rlat)

   do ii = 1,ncols
      i = ind(ii)
      rhocgs = rho(i)*1.e-3_r8     ! density in cgs units
!
! exponential temperature factor
!
!        efact = exp(0.025*(t(i,k)-t0))
!
! some temperature dependent constants
!
!!$      wt = min(1._r8,max(0._r8,(t0-t(i,k))*rdtice))
      wt = fice(i)
      icrit = icritc*wt + icritw*(1-wt)
!
! jrm Reworked droplet number concentration algorithm
      ! Start with pressure-dependent value appropriate for continental air
      ! Note: reltab has a temperature dependence here
      capn = capnw + (capnc-capnw) * min(1._r8,max(0._r8,1.0_r8-(p(i,k)-0.8_r8*p(i,pver))/(0.2_r8*p(i,pver))))
      ! Modify for snow depth over land
      capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,snowh(i)*10._r8))
      ! Ramp between polluted value over land to clean value over ocean.
      capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,1.0_r8-landm(i)))
      ! Ramp between the resultant value and a sea ice value in the presence of ice.
      capn = capn + (capnsi-capn) * min(1.0_r8,max(0.0_r8,seaicef(i)))
! end jrm
!      
#ifdef DEBUG2
      if ( (lat(i) == latlook(1)) .or. (lat(i) == latlook(2)) ) then
         if (i == ilook(1)) then
            write(iulog,*) ' findmcnew: lat, k, seaicef, landm, wp, capnoice, capn ', &
                 lat(i), k, seaicef(i), landm(i), wp, capnoice, capn
         endif
      endif
#endif

!
! useful terms in following calculations
!
      rat1 = rhocgs/rhow
      rat2 = liqmr(i)/capn
      con2 = (rat1*rat2)**0.333_r8
!
! volume radius
!
!        r3l = (rhocgs*liqmr(i)/(1.333*pi*capn*rhow))**0.333 * 0.01 ! meters
      r3l = con1*con2
!
! critical threshold for autoconversion if modified for mixed phase
! clouds to mimic a bergeron findeisen process
! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i)))
!
! autoconversion of liquid
!
!        cwaut = 2.e-4
!        cwaut = 1.e-3
!        lcrit = 2.e-4
!        lcrit = 5.e-4
!        pwaut = max(0._r8,liqmr(i)-lcrit)*cwaut
!
! pwaut is following tripoli and cotton (and many others)
! we reduce the autoconversion below critpr, because these are regions where
! the drop size distribution is likely to imply much smaller collector drops than
! those relevant for a cloud distribution corresponding to the value of effc = 0.55
! suggested by cotton (see austin 1995 JAS, baker 1993)

! easy to follow form
!        pwaut = capc*liqmr(i)**2*rhocgs/rhow
!    $           *(liqmr(i)*rhocgs/(rhow*capn))**(.333)
!    $           *heavy(r3l,r3lcrit)
!    $           *max(0.10_r8,min(1._r8,prlloc(i)/critpr))
! somewhat faster form
#define HEAVYNEW
#ifdef HEAVYNEW
!#ifdef PERGRO
      pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * &
              max(0.10_r8,min(1._r8,prlloc(i)/critpr))
#else
      pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* &
              max(0.10_r8,min(1._r8,prlloc(i)/critpr))
#endif
!
! autoconversion of ice
!
!        ciaut = ciautb*efact
      ciaut = ciautb
!        psaut = capc*totmr(i)**2*rhocgs/rhoi
!     $           *(totmr(i)*rhocgs/(rhoi*capn))**(.333)
!
! autoconversion of ice condensate
!
#ifdef PERGRO
      psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut
#else
      psaut = max(0._r8,icemr(i)-icrit)*ciaut
#endif
!
! collection of liquid by rain
!
!        pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994)
      pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton)

      pracwo(i)=pracw

!!      pracw = 0.
!
! the following lines calculate the slope parameter and snow mixing ratio
! from the precip rate using the equations found in lin et al 83
! in the most natural form, but it is expensive, so after some tedious
! algebraic manipulation you can use the cheaper form found below
!            vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs)
!     $               *0.01   ! convert from cm/s to m/s
!            snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i))
!            snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04
!            lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25
!            csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd)
!
! coefficient for collection by snow independent of phase
!
      csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05

!
! collection of liquid by snow (lin et al 1983)
!
      psacw = csacx*liqmr(i)*esw
#ifdef PERGRO
! this is necessary for pergro
      psacw = 0._r8
#endif

      psacwo(i)=psacw

!
! collection of ice by snow (lin et al 1983)
!
      psaci = csacx*icemr(i)*esi
!
      psacio(i)=psaci

! total conversion of condensate to precipitate
!
      ptot = pwaut + psaut + pracw + psacw + psaci
!
! the recipricol of cloud water amnt (or zero if no cloud water)
!
!         rcwm =  totmr(i)/(max(totmr(i),small)**2)
!
! turn the tendency back into a loss rate (1/seconds)
!
      if (totmr(i) > 0._r8) then
         coef(i) = ptot/totmr(i)
      else
         coef(i) = 0._r8
      endif

      if (ptot.gt.0._r8) then
         fwaut(i) = pwaut/ptot
         fsaut(i) = psaut/ptot
         fracw(i) = pracw/ptot
         fsacw(i) = psacw/ptot
         fsaci(i) = psaci/ptot
      else
         fwaut(i) = 0._r8
         fsaut(i) = 0._r8
         fracw(i) = 0._r8
         fsacw(i) = 0._r8
         fsaci(i) = 0._r8
      endif

      ftot = fwaut(i)+fsaut(i)+fracw(i)+fsacw(i)+fsaci(i)
!      if (abs(ftot-1._r8).gt.1.e-14_r8.and.ftot.ne.0._r8) then
!         write(iulog,*) ' something is wrong in findmcnew ', ftot, &
!              fwaut(i),fsaut(i),fracw(i),fsacw(i),fsaci(i)
!         write(iulog,*) ' unscaled ', ptot, &
!              pwaut,psaut,pracw,psacw,psaci
!         write(iulog,*) ' totmr, liqmr, icemr ', totmr(i), liqmr(i), icemr(i)
!         call endrun()
!      endif
   end do
#ifdef DEBUG
   i = icollook(nlook)
   if (lchnk == lchnklook(nlook) ) then
      write(iulog,*)
      write(iulog,*) '------', k, i, lchnk
      write(iulog,*) ' liqmr, rainmr,precab ', liqmr(i), rainmr(i), precab(i)*8.64e4_r8
      write(iulog,*) ' frac: waut,saut,racw,sacw,saci ', &
           fwaut(i), fsaut(i), fracw(i), fsacw(i), fsaci(i)
   endif
#endif

   return
end subroutine findmcnew

!-----------------------------------------------------------------------------
! Sets rhu to a different value poleward of +/- 50 deg latitude and 
! levels above the tropopause if cldwat_polstrat_rhmin is specified
! ** This is used only for special waccm/cam-chem cases with cam4 physics **
!-----------------------------------------------------------------------------
subroutine relhum_min_adj( ncol, troplev, dlat, rhu, rhu_adj )

  integer,  intent(in)  :: ncol
  integer,  intent(in)  :: troplev(:)
  real(r8), intent(in)  :: dlat(:)        ! latitudes in degrees
  real(r8), intent(in)  :: rhu(:,:)
  real(r8), intent(out) :: rhu_adj(:,:)

  integer :: i,k

  rhu_adj(:,:) = rhu(:,:)
  if ( .not.do_psrhmin ) return

  do k = 1,pver
     do i = 1,ncol
        if ((k .lt. troplev(i)) .and. &
             ( abs( dlat(i) ) .gt. 50._r8 ) ) then
           rhu_adj(i,k) = psrhmin
        endif
     enddo
  enddo

end subroutine relhum_min_adj

end module cldwat
